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1. Introduction and overview 

The inspiraling of black hole binaries is receiving much recent attention both because it is 
an exciting potential source of detectable gravitational waves|]IJ], and due to its inherent 
interest as a strong field gravitational interaction. The process of inspiral and merger 
is being investigated with a number of techniques. Newtonian and post-Newtonian 
computations are appropriate to the early stages of inspiral; numerical relativity || 
and black hole perturbation theory[|J are used for the late strong field stage of the 
inspiral. The present paper deals with aspects of the intermediate phase, when strong 
field effects are too important for post-Newtonian methods to be useful, but in which the 
full power of numerical relativity is not required. Detweiler and collaborators have 
drawn attention to the possible value of an approximation based on the comparison 
of the orbital time r or b for the binary holes with the time r ra d on which gravitational 
radiation acts to change the orbit. From a dimensional analysis of an equal mass binary 
of mass M and separation a this ratio is 

z=i oc (^KT . (i) 

Trad \ CL J 

The factor GMj (ac 2 ) on the right is an indicator of "how relativistic" the gravitational 
interaction is between the two holes. When a is on the order of 30 or so times GM/c 2 , 
the ratio in ([I]) will be small, and the orbits may be quasiperiodic, but nonradiative 
relativistic effects may still be important. 

It is useful to focus on a particular strong field relativistic effect of great importance: 
the innermost stable circular orbit (ISCO). The particle limit in gravitation theory treats 
the mass /i of one of the orbiting objects as much smaller than M, the mass of the other. 
In this approximation, there is a minimum radius for circular orbital motion, at a value 
r = 6GM/c 2 in the case of a nonrotating hole. The existence of this limiting orbit is a 
purely relativistic effect; there is no such limit in Newtonian theory. It is, furthermore, 
unrelated to radiation. In the particle limit radiation reaction is a force of order (/x/M) 2 
and is ignored; the particle moves on a geodesic. There is no such justification for 
neglecting radiation reaction for the case of binary motion of equal mass holes. Since all 
orbits are being degraded by gravitational radiation, there is no meaning to "stability" 
in principle. But there is an important practical question: Do relativistic forces arise 
in the late binary motion that drive orbiting particles to plunge inward on a time scale 
much shorter than the time scale due to radiation? It is very plausible that this sort 
of "practical" meaning can be given to the question of stability. For radial infall of 
equal mass holes, we know|| that the radiated energy is a very small fraction of the 
mass energy of the system as the holes fall into each other from moderate separations. 
Gravitational radiation reaction, therefore, cannot be a significant modification of the 
motion of the holes. 

This and other questions can be investigated in the absence of radiation reaction 
by seeking an approximation to the slowly evolving spacetime which is periodic, or, 
in the case of circular orbits, stationary. It turns out that this greatly changes the 
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mathematical nature of the solution process. The problem of evolving Cauchy data is 
converted into a boundary value problem. Past experience with this type of problem 
indicates that this boundary value problem suffers from none of the instability difficulties 
of numerical nonlinear evolution. The investigation of such questions has usually been 
based on an ansate|7|, |8| for suppressing radiative degrees of freedom in general relativity 
(GR), but such approaches by necessity only solve a subset of the full Einstein equations. 
We propose a very different approach. We do not attempt to suppress radiation fields, 
rather we suppress radiation reaction by requiring that the energy lost to outgoing 
radiation be replaced by a corresponding amount of ingoing radiation. 

Since the goals of this paper are to introduce the general ideas behind a new 
approximation scheme, and to demonstrate the numerical implementation of this 
scheme, we devote our attention here to a toy model rather than to GR, with its added 
complexities. The toy model is the simplest system containing the relevant features of 
radiation, nonlinearity, and decaying orbits, namely that of a nonlinear scalar field in 
2+1 dimensions. The choice of two rather than three spatial dimensions is made for two 
reasons: First and foremost, it makes the problem of finding a stationary solution less 
computationally intensive, since the wave equation has to be solved on a two- rather than 
three-dimensional grid. Second, a scalar field theory in 2+1 dimensions is equivalent 
to the same theory in 3+1 dimensions with all the sources and fields required to be 
translationally invariant in one spatial dimension. This theory with line-like sources is 
analogous to the problem of orbiting line-like sources in 3+1 GR. (Note that the scalar 
theory in 2+1 dimensions is not analogous to a problem in 2+1 GR. In 2+1 GR, gravity, 
i.e., Riemann curvature, vanishes outside gravitating bodies.) The formalism for 3+1 
line-like sources in GR has already been developed [[H]], and the spacetime they generate 
will be the topic of a subsequent paper. 

The analogy, with either point-like or line-like sources, between the scalar field 
problem and full GR is of course only a qualitative one, and we consider the scalar 
field results only as a test of the general method and not in any way a realistic model 
of the problem in GR. By the same token, the study now in progress of orbiting line 
sources in GR will not serve as a realistic simulation of the astrophysical problem with 
localized sources, but rather as a toy model that incorporates the added complexity 
of a gravitational problem while remaining computationally similar. Line-like sources 
in GR have features not found with localized sources, such as the lack of an ISCO 
for test particles and the lack of an asymptotically flat region in which gravitational 
waves have familiar propertiesQ. On the other hand, the analogy between the 2+1 
and 3+1 problems in scalar field theory is much closer, as illustrated analytically in 
Appendix B| in the absence of nonlinearities. We therefore consider the difference 



between looking numerically at 2+1 rather than 3+1 scalar field theory to be one 
primarily of computational complexity. Finally, 

Our toy model consists of two particles, each with a charge that couples to the 
nonlinear scalar field. To understand our method it is useful to consider three different 
solutions for the particle motion and the scalar field: I. The radiation is outgoing, and 
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as a result of the loss of energy due to the radiation, the orbiting particles spiral inward. 
II. The radiation is outgoing, but due to constraining forces the particles remain in 
circular periodic orbits. III. Again, the particles move in circular periodic orbits, but 
now the scalar radiation is balanced and the waves are some generalization of standing 
waves; there is as much ingoing radiation energy flux as outgoing. The solution of type 
I is the (scalar field, 2+1) analog of the problem of binary inspiral in GR. The solution 
of type II is a reasonable approximation of the type-I problem when T orb /r rad <C 1. 
This type of solution makes sense for a scalar field model; the constraining forces that 
maintain the periodic motion can be invoked ad hoc. Such forces need not couple to 
the scalar field, so they can be specified to have no effect on the problem other than to 
maintain the periodic circular orbits. In GR, on the other hand, all interactions couple 
to gravitation, and a solution of type II does not make sense. Solutions with periodic 
orbits and outgoing radiation should therefore not exist in GR. Solutions of type III, 
however, are not a priori ruled out by such considerations. 

There are two principal goals of this paper. The first is to show that the type-II 
fields can be found by solving a boundary value problem, and to suggest that such a 
problem is much more easily solved than a Cauchy evolution problem. Our second goal 
is to demonstrate that a solution of the type-Ill problem can also be computed without 
evolution, and that the computed solution gives a good approximation to the fields of the 
type-II problem, when the conditions of our approximation are valid. The implication 
is that the physical problem in GR, the type-I problem, can be approximated from the 
"easily" solved type-Ill problem. 

The remainder of this paper will be organized as follows. The basic toy model, a 
nonlinear scalar field will be introduced in Section |2] and solutions, both analytic and 
numerical, will be presented for periodic fields with outgoing boundary conditions. In 
Section ^| we discuss the meaning of radiation balanced fields and we make a particular 
choice (the "TSGF") of such fields for our nonlinear model. Issues related to the 
application of the quasi-stationary approximation to GR are discussed in Section |. 
Appendix A gives some technical details of the general solution to the linear wave 
equation with equal flux of ingoing and outgoing radiation. |Appendix B| presents the 
solution to a linear scalar field problem in 3+1 dimensions, as an illustration of the 
relationship of the 3+1 and 2+1 problems. Finally, Appendix C demonstrates explicitly 
(in the linear case) the vanishing of the radiation reaction force for the radiation- 
balanced solutions. 



2. Nonlinear scalar field with periodic orbits 

Our toy model is based on a nonlinear field ip in 2+1 dimensions satisfying 

V 2 ^ - dfy + AJF(V>, p) = -a . (2) 

where a is a scalar source for the field. The spacetime for the field is Minkowskian 
and we use polar spatial coordinates p, <fi. The term \T(il),p) is a term nonlinear in ip 
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that may be an explicit function of p, but not an explicit function of or t. For later 
computational convenience we require that the nonlinear term satisfy the symmetry 
condition 

;F(V,p) = -.FHM • (3) 

The constant A is a parameter that governs the strength of the nonlinear term. We now 
specialize to the source 

a = a - 1 Q5(p - a) [8((f> -Qt- tt/2) - <5(0 - M - 3rr/2)] (4) 

representing two point sources (in two spatial dimensions) of scalar charge density ±Q, 
moving around each other in circular orbits of radius a, at angular speed Q. 

If we were to treat (J|) as an evolution problem, we would have to specify Cauchy 
data and integrate forward in time to find ip(p, 0, t). Instead we concern ourselves only 
with the steady state solution, a solution that would evolve after transients associated 
with the initial conditions are radiated away. This steady state solution would embody 
the symmetries of the source. In particular, the source depends on and t only in the 
combination ip = — Qt and we seek a solution with the same symmetry. That is, we 
look for a solution ip(p,(p). Since the left hand side of (0) does not have an explicit 
dependence on or t, such a solution is allowed. 

When (||) is restricted to solutions of the form ip(p, ip) it reduces to 



ld_ ( df\ 
pdp V dp J 



- 2 -& 
P 

a- l Q8(p 



d 2 i) 

dip 2 



+ XFty,p) = -a(p,(p) 



a) [% - 3vr/2) - % - tt/2)] 



(5) 
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Figure 1. Coordinate regions for the periodic wave equation. The sources are shown 
as dark circles. The light cylinder is shown as the line with dots and dashes at p = 1/57. 
The outer boundary, p = /0 max , is not shown. 

The coordinate regions for this equation are shown in Figure |I[ Although this 
figure misrepresents the topology (a disk) of the physical problem, it is the appropriate 
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description of the (p, ip) grid on which ([|) is to be solved as a finite difference equation. 
The specific conditions to be used in solving the problem are the following: (i) Due to 
the symmetry of our source, and to condition (|3|), the solution is to have the symmetry 
ip(p,p) = —tp(p,p + 7r). This allows us to restrict the numerical solution to the range 
< (p < 7r. (ii) Due to the antisymmetry under <p — > <p + 7T, the solution for ip must 
vanish at all p — grid points, (iii) Outgoing Sommerfeld boundary conditionsare 
imposed at some maximum radius p max by requiring 

/ dip 



dp dtp 



= . (6) 

p— pmax 

If (|5|) and the above conditions (i-iv) are put on a grid of N r radial lines and N v lines, 
containing N = N r x N v grid points with a priori unknown values of ip, a system of 
equations for these unknowns is found. The solution of this system is the finite difference 
solution of our physical problem. 

The finite difference procedure just outlined is relatively straightforward, but it has 
an unusual feature. Note that the nature of the differential equation in (gj) changes 
at the "light cylinder" p = I/O, shown as a line with dots and dashes in Figure |]. 
For p < I/O the equation is formally elliptic, while for p > I/O it is hyperbolic. 
Typically elliptic partial differential equations are solved as boundary value problems, 
with auxiliary data given on a closed boundary surrounding the region of the solution, 
but hyperbolic equations are given Cauchy data on an "initial" hypersurface. The 
common wisdom is that a hyperbolic problem with data specified on a closed surface 



can have more than a single solution p2| . Despite this we treat the entire coordinate 
region (0 < p < tc, < r < p max ) as a boundary value problem. We have not attempted 
to give a rigorous proof that the boundary value approach to @ and (|J) is well posed, 
but several nonrigorous justifications are worth mentioning: (i) Although nonuniqueness 
is a possible feature for a hyperbolic equations with boundary values, whether or not 
a particular problem suffers from this difficulty depends on details of the problem, not 
only on whether it is hyperbolic, (ii) The physical problem described by the boundary 
value problem appears to be well posed, (iii) Numerical solutions of the boundary value 
problem are stable and insensitive to numerical grid size. 

Here, as in the next section, it will be useful to separate the complexities of 
nonlinearity from other issues. To do this we temporarily set A to zero so that (|) 
becomes a linear equation. The outgoing radiation solution for this problem is easily 
found, with standard techniques, in the form of a series of Bessel J m and Neumann N m 
functions. For p > a, this series is 

Vw = Q X! (— l)( m+1 )/ 2 J m (mVla) [N m {mVtp) sin mp — J m {mVtp) cos mp] . (7) 

m=l,3,5,... 

This solution shows, by example, that there are no hidden difficulties in finding a 
solution to (|5|) in the case A = 0, and hence there is no fundamental problem in 
using a boundary value approach to solve a problem with outgoing radiation. It also 
demonstrates explicitly that, at least in the A = case, the light cylinder p = I/O is 
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not a special surface in the problem. Note that ([/]) gives the solution for waves that 
are "outgoing at infinity." Solving the linear problem for waves that are "outgoing" at 
a finite radius, i.e. , for the boundary conditions in (||), is almost as simple as for true 
outgoing waves. For p > a, the solution is 

^ = T/Wt-Q £ (-l) (m+1)/2 J m (mfia)J m (mfip)Re( 7m e^) , (8) 



where 



In 



m=l,3,5,. 



Hg){z)+ifHg){z) 



Here H$ indicates the Hankel function of type 1. From a well known property of Hankel 
functions]!^] the numerator of @ vanishes as mf2p max — > oo. This demonstrates that 
(aside from roundoff and truncation error) a solution on a grid of finite radial extent 
approaches the true outgoing solution as the radial extent of the grid becomes infinite. 
We have also checked that there are no difficulties, with truncation or otherwise, in 
solving (|5|) with finite difference methods. The solution found in this way was compared 
with ([71). The agreement was excellent, and the expected second order convergence of 
the finite difference scheme was confirmed. 

To investigate a nonlinear model numerically a specific choice must be made for 
the nonlinear term T . To suit our purposes this term must satisfy several criteria in 
addition to the symmetry condition in (|3|). One criterion is that this nonlinear term be 
small enough at the outer boundary p max , so that the outgoing condition (||) is a good 
approximation at a large finite radius. The nature of this requirement can be seen if 
we choose T simply to be if), so that the nonlinear term in @ becomes \ip, and (||) is 
the Helmholtz equation. The "outgoing" solution to this linear problem is that given in 
(^), except that the arguments of the Bessel and Neumann functions have y/m 2 Q 2 + A 
in place of mQ. At large p this solution will satisfy the condition (|6]), for a particular 
angular Fourier mode, if the f2 is replaced by \/m 2 VL 2 + A jm. The standard outgoing 
condition is then a good approximation only if |A| <C Q 2 . To get a rough idea of the 
effect of nonlinearities on boundary conditions we can view the nonlinear term AJF as 
\ e f{ip with A e ff, the effective A, taken to be A(JF)/ (ip). Here "()" indicates some sort of 
average (perhaps an r.m.s average over all <p and one wavelength). A rough criterion 
for a nonlinear term that is compatible with the boundary condition (|6]) is that A c g, at 
Pmax, be much smaller than Q 2 . For our toy model to have interesting nonlinear effects, 
however, there is a somewhat contradictory requirement: the nonlinear term must be 
significant, even dominant, at small radius. A nonlinear term like T = ip 3 can be strong 
at small radii and weak at large radii if the field ip falls off quickly enough. For our line 
sources, the radiation fields fall off only as r _1//2 . This slow fall off causes computational 
difficulties in the application of boundary conditions. For this reason we include a 
factor exp (— (ap/a) 2 ) in the nonlinear term. To be certain that the nonlinear source is 
sufficiently well behaved near the point sources at the points (a, vr/2) and (a, 37r/2), we 
choose a source term that does not diverge at those points. The specific choice used in 
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our numerical investigations has been 

T = exp(-(ap/a) 2 ) — — . (10) 
(ip 2 + 1) ' 

The outgoing wave pattern for ip is illustrated in Figure || The field ip is shown 
as a function of pcos0 and psin0 for two different times: (a) fit = 0, and (b) fit = n. 
The figure illustrates how the rotation of the "rigid" pattern sends waves outward. The 
results were obtained from numerical runs with 961 x 41 gridpoints in p G [0,40] and 
<p G [0, 7r]. The parameter values are Q — 1, ft — 0.5, a = 0.5, A = 1500, a = 1. At 
fit = 0, the positively charged particle is at = 7r/2 and at Qt = it it is at <fi = 3tt/2. 
Figure [3| shows the importance of nonlinear effects for the outgoing fields of Figure 0. 





Figure 2. The field ip as a function of pcos<p and psmcj) for the two different 
times fit ~ (a) and Qt = n (b). This is a numerical solution to the nonlinear wave 
equation (with nonlincarity parameter A = 1500 and localization parameter a = 1) in 
the presence of outgoing-radiation boundary conditions. The sources are located at a 
radius of p = 0.5. Note that there are no discernable irregularities at the light cylinder 
p = 2. 



The field ip is shown as a function of p at (p — 7r/2, for Q = 1, fl = 0.5, and a = 0.5. 
The solid curve shows the solution of the linear problem for A set to zero; the dashed 
curve shows the nonlinear fields for A = 1500, a = 1. Details of the fields near the 
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source particles are shown in the insert. Though the nonlinear terms are negligible in 
the wave zone, the amplitude of the waves is more than doubled by the increase in the 
effective source strength produced by the nonlinearity in the central region. 




Figure 3. The linear (solid curve) outgoing solution for A = and the nonlinear 
(dashed curve) outgoing solution for A = 1500, a = 1. Both curves are for Q = 1. 
£1 = 0.5, and a — 0.5. The fields are shown at p — ir/2. The sharp peaks correspond to 
the location of the particle, at p = 0.5. Note that there are no discernable irregularities 
at the light cylinder p = 2. 



3. Radiation-balanced boundary conditions 

The universality of gravitational coupling in GR means that periodic orbits and outgoing 
radiation are not compatible. In this section we investigate a choice of solution for 
our toy model that does seem applicable to periodic orbits in GR. We introduce 
boundary conditions that do not carry net energy away from the inner region of the 
orbiting objects, boundary conditions containing equal measures of ingoing and outgoing 
radiation. There is a temptation to relate this idea to that of "standing waves," and 
to impose Dirichlet boundary conditions that if) vanishes at some "wall" located at 
p = Pwaii- As in the previous section, it is useful to separate issues of nonlinearity from 
other issues, and to look at the A = case first. For this linear problem the solution, 
for p > a, with ip vanishing at p = p wa u is 

4> = Q E (-l) (m+1)/ V m (mfia) [N m (mQp) + P m J m (mQp)) sin rwp , (11) 

m=l,3,5,... 

where 

iV m (mQ/wi) 

Pm — t / o \ • 

J m (miZp wa ii) 

The expressions in flTTp and ( |T^ ) do not describe a well behaved solution. In general, 
the denominator for (3 m will come arbitrarily close to zero, as larger and larger values 
of m are included in the sum. We have confirmed that the finite difference solution to 
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(I!]) and fll2|) is not stable for small changes in the location of p wa ii, or for a change in 
the number of grid points. [Note that increasing the number of angular grid divisions is 
roughly equivalent to increasing the maximum value of m included in the sum in ((TTJ) .] 
In considering fields ip that are acceptable mixtures of ingoing and outgoing waves, 
it is again helpful to look at the linear (A = 0) problem. In this case, a more-or-less 
obvious acceptable choice of ip can be constructed simply by averaging the solutions to 
(|) with ingoing waves at infinity and with outgoing waves at infinity. The result, for 
p > a, is 

tlj = Q (-l) {m+l)/2 J m (mtta)N m (mVLp)sinmp. (13) 

m=l,3,5,... 

Although this is not the most general solution of the linear problem with equal amounts 
of in and outgoing waves, it is the most natural solution, as discussed in [Appendix A[ 

The field in ([13|) is a solution to @ in the A = case, but not for any simply 
stated boundary conditions. In particular, it does not correspond to the vanishing of 
ip at some specific finite radius. Nor is it the limit in which the "wall" of (|TTD, ( |T2| ) is 
at oo; no such limit exists. Rather, it is a superposition of two solutions each of which 
has a simply stated boundary condition (ingoing or outgoing), and superposition is not 
valid for solutions of nonlinear equations. We therefore reformulate our toy model so 
that fields can be found that are the nonlinear equivalent of (|l~3|) . In place of a partial 
differential equation, we introduce an integral equation. The first step in doing this is 
to rewrite (|) as 

Lip = a eS (p,p,ip) , (14) 



where 



and 



c = --— A _ (1 _ n ^ _*!_ (is) 

pdp dp Ip 2 I dp 2 



a eff (p, <p, ip) = a(p, <p) + \T{ip, p) . (16) 
We take the solution of ( |l4"[jT6|) to be the time symmetric Green function (TSGF) 



solution given by 

^tsgf(p,<p) = / / G T s(p,V,P,v')^es(p ,<p'ip(p',(p'))p'dp'd<p' , (17) 



where the Green function Gts(p ; V 9 ! p' i V 9 ') is the time symmetric inverse to the linear 
operator £ corresponding to equal mixtures of ingoing and outgoing waves. In principle, 
this Green function can be written explicitly as 



G(p, <p; p', <p')ts 



1 1 00 

-— ln(p/p') - - V J m (mttp')N m (mQp) cosm(<p - ip') p> p' 

2n 2 m=l 

oo 

- N m (mVLp')J m (mVtp) cosm(<p - ip') p < p' 

* m=l 

(18 



Quasi- stationary binary inspiral II: Radiation-balanced boundary conditions 



11 



In practice, this series form of the time symmetric Green function is not directly 
applicable to our numerical method, and it is more useful to write the same radiation 
balanced solution with the following symbolic notation. We denote the solution to the 
nonlinear problem with outgoing boundary conditions as 

^out = £out°"eff , (19) 

where £~ u \ is the inverse to C for outgoing boundary conditions, i.e. , the retarded- 
time Green function. The ingoing solution ip- m is written in a parallel manner using 
the advanced-time Green function C\ n . The linear superposition of the ingoing and 
outgoing (LSIO) solutions is not itself a solution since the problem is nonlinear. 

To arrive at a field ■j/'tsgf that is a solution, and that corresponds to the solution in 
(|17|), we superpose the operators £j~ and £~ ut and take our radiation balanced solution 
to be 

^TSGF = g (Cfr + £out) a cS = ^TSGF^eff- (20) 

This integral equation is to be solved by iteration. From t/'tsgf,™, the n th iteration for 
^tsgf, we construct cr e s >n , and then find the (n + l) th iteration of the solution from 

V'TSGF^+l = ^TSGF^cfT.n • (21) 

This method is well suited to implementation as a finite difference solution to a boundary 
value problem on a grid of iV points, like that in Figure [I]. In this implementation, t/'tsgf 
and cr e ff are one dimensional vectors of length N, and £, along with the chosen boundary 
conditions, forms a.n N x N matrix. Since the form of this matrix depends on boundary 
conditions, the numerical inverse C~ l is also specific to the boundary conditions. The 
numerical radiation balanced operator £tsgf i s simply the average of the matrix inverses 
of C~ l for the ingoing and outgoing problems. Results of this numerical solution are 
shown in Figure This figure, when compared with Figure [2] nicely illustrates the 
symmetry of the radiation field with respect to reversal of <p. It is not so effective in 
illustrating the fact that there is no radius at which ip vanishes at all values of (p. This 
is due to the strong dominance of the m = 1 multipole of the field. For larger values 
of the source velocity v = afl the m = 3, 5, . . . multipoles make stronger contributions, 
but for those values of v for which we could get accurate solutions (up to v around 0.8) 
the m = 1 multipole continued to dominate the appearance of the fields in a plot like 
that in Figure f|. 

We now argue that within the scope of our approximation, the TSGF solution in the 
wave zone is close to a linear superposition of the ingoing and outgoing (LSIO) solutions 
for the same orbiting sources. It is important to recall that due to the nonlinearity, 
the LSIO is not a solution of @, but it does have the convenient property that if 
the nonlinearities are weak in the wave zone, the outgoing and ingoing solutions can 
easily be extracted from the LSIO. We are claiming then that an approximate outgoing 
solution to the problem, in the wave zone, can be found from the TSGF solution. We 
emphasize that we are not claiming that nonlinear effects are weak. We will show that 
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Figure 4. The field i/j as a function of pcos</> and ps'mcj) (at time t — 0). This is 
the numerical solution to the nonlinear wave equation (with nonlincarity parameter 
A = 1500 and localization parameter a — 1) according to the TSGF method of 
imposing a radiation-balanced boundary condition. The sources are located at a radius 
of p = 0.5. Note that there are no discernible irregularities at the light cylinder p = 2. 



this approximation method is successful for problems with very strong nonlinear effects. 
The reason for this success can be seen most easily if fl2~0|) is rewritten as 



The TSGF solution and the LSIO superposition differ due to the nonlinear terms 
contained within the effective source. Those nonlinear terms will be significant only 
in the small radius inner regions of the physical space, where the fields are strong. But 
in the inner, strong field, regions the solution should not be highly sensitive to whether 
ingoing or outgoing boundary conditions are imposed. Indeed, this last statement is a 
way of viewing the underlying idea in our approach. If the fields near the orbiting objects 
are not significantly influenced by the distant boundary conditions, then radiation 
reaction cannot be important. We are now assuming the converse: with parameters 
for which radiation reaction is not important, the fields near the sources will not be 
sensitive to the boundary conditions. 

If the nonlinear contributions to cr c ff are nonnegligible only in the strong field region 
near the orbiting objects, and if the fields there are insensitive to boundary conditions, 
then we can conclude that o- cS (ip in ) and cr eff (-?/> out ) differ negligibly. It is then plausible 
that they are also negligibly different from cr e ff (^tsgf), and hence that the LSIO is 
approximately the same as the TSGF solution. This conclusion, furthermore, should be 
valid to the same degree that it is valid to ignore radiation reaction (more specifically 



V'TSGF = -A n V cff (^TSGF) + ^ ^out^eff OtSGf) 

and compared with the LSIO (not a solution of the nonlinear problem), 



(22) 




(23) 
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to ignore the sensitivity of the fields near the source to the boundary conditions on the 
waves) . 



For the nonlinearity given by (10), with A = 1500, a = 1, it turns out that the 
approximation of TSGF by LSIO is too good for an effective illustration. In this case 
a plot shows no discernible difference between the TSGF and LSIO, even though the 
nonlinearity plays a strong role, as shown in Figure [| A discernible difference can 
be seen if the a parameter is reduced. The effect is to spread the region of strong 
nonlinearity to larger radius where the fields are somewhat sensitive to the boundary 
conditions. The comparison of the TSGF and the LSIO fields, for a = 0.1 is shown in 
Figure^. The field ip is shown as a function of p both for (p = ir/ '2, the angular position 
of a source particle, and at tp = tt/A. No difference between the TSGF and the LSIO 
fields can be seen at small radii. In the radiation zone a small phase shift can be seen 
between the LSIO and the TSGF and the LSIO waves are seen to have a slightly smaller 
amplitude. 

In our approximation approach, the idea is to find the solution of the outgoing wave 
amplitude from the TSGF solution, by considering the TSGF solution in the wave zone 
to be (approximately) an equal mixture of ingoing and outgoing radiation. The specific 
application of this idea requires fitting the TSGF radiation field to a sum of multipoles 
of the form 

^tsgf ~ C m N m (mQp) smrmp, (24) 

m=l,3,5... 

and from the amplitudes C m inferred from the fit, writing the outgoing solution as 

V'out ~ X! C m [N m (mQp) smmip — J m (mQp) cos mip] . (25) 

m=l,3,5... 

Using our toy nonlinear model, we can solve for both the TSGF and the outgoing 
solution, to check the accuracy of this method. For the model with parameters Q = 1, 
a = 0.5, Q = 0.5, A = 150, a = 0.1, fitting the TSGF from p = 20 to p = 40, gives a 
wave amplitude for outgoing radiation that is larger than the true amplitude by 24%. 
For a = 1 the discrepancy is only 0.04%. 



4. Conclusion and discussion 



For our nonlinear toy scalar field model, we have demonstrated that a radiation balanced 
field, the time symmetric Green function (TSGF) solution, can be found by solving a 
numerical problem similar to a boundary value problem. Although the character of 
the differential operator is elliptical inside the light cylinder and hyperbolic outside, 
no special treatment of this surface was necessary, and the solution was found with 
boundary value methods usually associated with elliptical equations, despite the outer 
boundary being in the "hyperbolic" region. We have also shown that to the extent that 
radiation reaction forces can be ignored, the TSGF solution is a good approximation to 
the solution with outgoing boundary conditions. 
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Figure 5. The difference between TSGF and LSIO. The figures show the source 
region and the wave zone for a = 0.5, il = 0.5, a = 0.1 and A = 150. The plot on the 
right is for ip = 7r/4 and that on the left is for ip — w/2, the ip location of one of the 
particles. 

We must now ask how related methods might be brought to bear on the problem 
of orbiting objects in GR. An important feature of the nonlinear scalar toy model is 
that the nonlinearities in the theory do not occur in the wave operator. This allowed us 
directly to recast the problem as an integral equation by using the time-symmetric Green 
function for the wave operator. It is not at all clear how the GR equations can be cast 
in a form with a linear wave operator, or whether it is in principle impossible. But it is 
not necessary that we follow these same steps in dealing with the GR equations. What 
is required, most generally, is any method for specifying a time symmetric solution. 

For such a method in GR it is expected that some features of the basic physics will 
be the same as in the scalar model. In particular, for motions r or b/T ra( i <C 1 the nonlinear 
terms in the effective source should be insensitive to boundary conditions. It can then 
be supposed that the radiation balanced solution is a good approximation to the linear 
superposition of ingoing and outgoing solutions (LSIO) and that from the LSIO we 
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can infer the outgoing solution. There is, however, an important difference between this 
claim as it applies to GR, and the (demonstrably true) claim for the scalar field model. In 
the scalar field model we were comparing LSIO fields for periodic motion with the TSGF 
solution for periodic motion. In GR, there can be no LSIO solution for periodic motion, 
since ingoing, or outgoing, energy flux would be incompatible with periodic motion. 
For GR, the analogy needs to be made directly to the solution with outgoing waves 
and a slow rate of orbital decay due to radiation reaction. By constructing a spacetime 
with only outgoing radiation from the TSGF solution, we arrive at an approximation in 
which the radiation fields are periodic, with a period that is constant, not slowly drifting 
in time. In the notation introduced in Section [l], the quasi-stationary method in GR 
requires that the type-Ill problem be solved and used to produce an approximation to 
the type-I physics. 

An additional complication is that the spacetime of the TSGF fields cannot be 
asymptotically flat in GR, since there is an infinite amount of energy contained in the 
radiation fields JTT]] . This apparent pathology can be viewed as irrelevant if we think 



of our TSGF, or outgoing, solutions as being an approximation only in a finite central 
region of the space. It is worth noting that this viewpoint is consistent with the nature 
of the numerical solution that is based on boundary conditions applied at finite radius. 
The important question that remains is whether the inner region, the region in which 
the problem is solved, is large enough to include a zone in which waves are weak and 
in which the ingoing and outgoing waves of the approximate LSIO can be disentangled. 
To investigate this, we note that the gravitational wave luminosity of the binary is of 
order L GW ~ (GM/ac 2 ) 5 (c 5 /G), where we are using the same notation as for ([[]). The 
gravitational wave energy Eqw contained within a sphere of radius p max is of order 
L GW p max /c, and hence the ratio of E GW to the orbital energy GM 2 /a of the binary, is 
of order 



E GW r fGM\ 3 p max 
GM 2 /a \ac 2 ) a ' 1 } 

The nature of our approximation requires that GMj (ac 2 ) be small, so we conclude that 
the gravitational wave energy contained within p max will be much less than the orbital 
energy, even for values of p max large compared to a. 

This last conclusion is important if we are to hope to use versions of our method to 
make inferences about the ISCO. By investigating the dependence of the binary energy 
on radius, we can study whether there is an instability like that for particles. The 
"energy" of the orbit must be computed as a surface integral at a large radius. If this 
integral were significantly influenced by the energy contained in the waves, conclusions 
about orbital stability would be suspect. 

Whether or not the methods of this paper can be applied to the ISCO question, the 
use of these methods in GR would provide an important addition to the tools needed 
to understand black hole inspiral. A direct and obvious use of the method would be to 
provide Cauchy data for numerical relativity. 
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Appendix A. General Radiation-Balanced Solution 



We present here some of the details of the general solution for the linear scalar 
field theory with equal flux of ingoing and outgoing radiation. The Green function 
G(p, p\ p', ip') for the linear problem satisfies 



1 d d 
p dp^ 'dp 



Q 2 



dp 2 



G(p,<p]p',<p') = -p 5{p - p')5{p - p') 



(A.1) 



When the p dependence is represented in a Fourier series, and the Green function is 
written as G(p, p\ p', <p') =Re[X) Gm(p] p', p')e tmip ] , the equation becomes 

1 



1 9 d 2 / o2 1 N 
pop dp \ p 2 . 



Q m (p; p',<p') 



2np 



-5{p - p')e 



■irrup 



(A.2) 



Reality of the Green function G(p,ip; p',<p') means that (?_ m (p; p', p') = G m (p; p' , ip')* , 
and the m = mode is irrelevant for the sources we consider, so we only need the 
general solution to ( A.2|) for m > 0, which is 



Q m (p;p',<p' 



--J m {mQp')N m {mQp)e- im ^' + T m {p', p')J m {mQp) p > p' 



N m (mnp')J m (mn P )e- im ^' + T m (j/, cp')J m (mQp) p < p' 



(A.3) 



For p > p', Q m can be written in terms of Hankel functions as 

G m = A m H^ (mfip) + B m HjV (mQp) , (A.4) 

where 

A m = ^r m (p',<p') - ±-j m ( m n P )e- im *' 

B m = \v m {p\ p>) + 1 J m (mfip)e- in ^ . (A.5) 

The condition for equal flux of ingoing and outgoing radiation is \A m \ = \B m \, for 
all m. This, and the expressions in ( |A.5| ), require that T m (p',ip') be a real multiple 
of J m (mVlp)e~ imLpl . The proportionality constant may be different for each m, so the 
general radiation balanced solution is 



g m (p;p',ip' 



--J m (mnp)N m (mn P )e- im ' p ' + g m J m {mnp')J m {mQp)e- im(p ' p > p' 
-]N m (mQp')J m (mnp)e- imip ' + g m J m (mnp')J m (mnp)e- imifi ' p < p' 
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Any choice of the set of real constants g m gives a radiation balanced solution. The 
choice made for the TSGF solution, g m = 0, for all m, is convenient for a numerical 
implementation that does not explicitly use Fourier decomposition. The choice g m = 
also seems to be the most natural one if one views the T m terms in ( |A.2| ) as free waves, 
disconnected from the source. 



Appendix B. Three-plus-One Dimensional Theory 

Here we present the solution of a linear toy model for point-like sources. We use standard 
spherical coordinates (r, 9, 0) and suppose that two particles of opposite scalar charge 
±Q orbit at frequency fi in the equatorial plane (6 = 7r/2), at radius a, with angular 
separation A0 = -k. In the linear differential equation of this toy model 

V 2 ^ - fi?V = -a~ 2 Q5(r - a)5(cos6) [5(0 - fit - tt/2) - 5(0 - fit - 3tt/2)] , (B.l) 

we make the ansatz that the solution is not of the general type ip(r, 8, 0, t), but rather 
of the type -0( r , 0, V 9 ) where = — fit. The differential equation then reduces to the 
point-like linear analog of (||): 

<9 2 ^ 



11^+— -fsin^U 
r 2 9r I <9r J r 2 sin 6 86 \ 89 J 



dip 2 



_r 2 sin 2 6 1 

= a- 2 Q5(p - a)5(cos 0) [8((p - tt/2) - 5(ip - 3vr/2)] . (B.2) 

If ip(r, 8, (p) is decomposed as a sum of spherical harmonics Ye m (6, <p) it is straightforward 
to find the solution that is well behaved at r = and that corresponds to outgoing waves 
at infinity. This solution, the point-like equivalent of (0), is 

^ f ji(m^la)hf (mfir) r > a 

iP(r,e,<p) = QnJ2 K t™Ye m (e,<p)x\ '/) (B.3) 

im yjt(m£lr)h\ '(mQa) r<a 



Here ji and lip indicate spherical Bessel and Hankel functions. The sum in (p.3|) is 
over odd i and odd m, and the coefficients K£ m are given by 

K lm = tm [F/ m (vr/2, tt/2) - F/Jtt/2, 3tt/2)] = (-l)( m - 1 )/ 2 2mF, m ( 7 r/2, 0) . (B.4) 
The point-like equivalent of (0), the radiation balanced TSGF solution of the problem, 



is 



[ jt{rn£la)ni(m£lr) r > a 
$(r,6 ) (p) = iQn^K tm Y tm (6,(p)x< . (B.5) 
« m [ je(mllr)ne(mlla) r<a 

where is the spherical Neumann function. These solutions serve as explicit 
illustrations that for point-like sources in 3+1 dimensions the field embodies the same 
symmetry as the sources. That is, the field rotates "rigidly." The dependence on time 
and on azimuthal angle appears only in the combination — fit. 
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Appendix C. Net Force on the Orbiting Particles 

In this appendix, we show that for a general radiation-balanced solution to the problem 
of orbiting particles, the proscribed circular orbits of the particles are consistent with 
the scalar field equations of motion without the need for external forces. (This is not 
true in the case of purely outgoing radiation.) We limit attention to the linear theory, 
so that we can easily separate out the forces due to each particle on the other, without 
worrying about any self-force. 

The covariant force law for a scalar-charged particle with charge Q, mass m, and 
energy-momentum vector p moving in a scalar field ip gives a force of 

m-yv/ = 4~ + m ~ lT uxP u P X = -QV^V • (CI) 
dr 

If we work in co-rotating coordinates, in which the particle is momentarily at rest, the 
metric is 

ds 2 = -(1 - Q 2 p 2 )dt 2 + dp 2 + p 2 d<p 2 + 2Qp 2 d<pdt (C.2) 

and the only non-zero component of the energy-momentum vector is 

p t = (1 - ftV)- 1/2 m =: 7 m , (C.3) 

so that the only relevant Christoffel symbol is T P t . (In the corresponding calculation in 
the 3+1 dimensional theory, only T r tt is relevant.) 

The equations of motion for p p and p^ in the 2+1 dimensional theory thus become 
dn p 

f = -Q^ + 7 2 mpfl 2 (C4) 
dr 

dp^ 

= -Qg™i/>„ , (C5) 

where we have used the form of the inverse metric in co-rotating coordinates and the 
fact that (di/j/dt) v = 0. The second term in ( C.4 ) represents the fictitious centrifugal 



force that the particle feels in the co-rotating reference frame. The radial momentum 
p p can remain zero if the repulsive effect of this term cancels the attraction due to the 
field, described in the first term. This implies a relationship among the orbital radius a, 
angular frequency Q, and charge-to- mass ratio Q/m. (The corresponding relationship 
is given for orbiting charged particles in the presence of a half-advanced, half-retarded 
electromagnetic potential in fI3| .) 



The angular force in ( |C.5| ), which represents radiation reaction force, will vanish if 
(and only if) the scalar field due to one particle has vanishing ip derivative at the location 
of the other particle. For concreteness, we consider the force on the positively-charged 
particle at p — a, cp — tt/2 due to the negatively- charged particle at p = a, = 3ir/2. 
The field due to the latter particle is 

if>(j>,<p) = -QG(p,<p;a,tor/2) , (C.6) 

where G(p, p>\ p', if') is the Green function used to construct the solution. Using the 
Fourier expansion of the Green function given in [Appendix A| , we find that 



ijj^(a, tt/2) = -Re 



(-J m {mtta)N m {mtta){-l) m + r m {a,3Tr/2)e im7r/2 



(C.7) 
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This vanishes if the quantity in large parentheses on the right-hand side is real; 
for a general radiation-balanced solution this is the case, since T m (a, 37r/2) = 
g m J m (mfla) exp(—im3ir / 2) with g m real. Thus the radiation-balanced solutions have 
zero radiation reaction force. 

However, in the case of purely outgoing radiation, T m (p',<p') = J m e~ imip ' /4z, and 
we can explicitly find the radiation reaction force as 

771 

Vv(a, tt/2) = - ]T(-l) m -[J m (mfta)] 2 ^ . (C.8) 

m 

A similar demonstration can be made in the 3+1- dimensional case using the 
outgoing-radiation and radiation-balanced fields ( |B.3| ) and ( |B.5| ), using the symmetries 
of the spherical Bessel functions under changes of sign in their arguments to illustrate 
that terms in the radiation-reaction force due to terms with opposite signs of m cancel 
each other out in the radiation-balanced mode expansion but not in the outgoing- 
radiation one. 
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